rng(5000, 'twister');

supply_est_yr = 2016; % the year in which we estimate entry costs

N_sim = 30;     % number of simulations


cutoff_month_count = 1; % count a brand as in a market if the brand is in a market for more than cutoff_month_count months

ind_yr = ismember(yr, supply_est_yr);

    
% use the same brewer_id and brewer_name for the merged brewers
ind = ismember(brewer_id, merge_id);
brewer_name_merged = unique(brewer_name(brewer_id == merge_id(1)));

brewer_id(ind) = merge_id(1);

brewer_name(ind) = brewer_name_merged;

brewer_loc_tbl.brewer(ismember(brewer_loc_tbl.brewer_id, merge_id)) = brewer_name_merged;

%% construct the set of potential products
% list of all brands
list_brand = unique(brand_descr(ind_yr));
n_brand = length(list_brand);

% brands that are not fixed
ind = craftdummy | ismember(brewer_name, {'shock top brewing co', 'michelob brewing co', 'blue moon brewing co'});

% brands that are fixed 
    % craft or macro beer firms' craft-like breweries
fixed_brand = setdiff(list_brand, brand_descr(ind))'; clear ind;

fixed_craft_brands = [];

for j = 1 : length(list_brand) % brands that are in less than 36 county/yearmonths
    if nnz(ismember(brand_descr, list_brand{j}))<=36
        fixed_brand = [fixed_brand, list_brand(j)];
    end
end



% get coordinates of all counties 
ind = ismember(geo_tbl.store_county_id, store_county_id(ind_yr));
list_county = geo_tbl.store_county_id(ind);
list_county_crdnt = [geo_tbl.latitude_cty(ind), geo_tbl.longtitude_cty(ind)];
n_county_cty = nnz(ind);

% compute attributes for products present in supply_est_yr
prod_attributes = cell(length(list_brand), 1);
is_fixed_prod = ismember(list_brand, fixed_brand);

for nb = 1 : length(list_brand)
    ind_nb = ismember(brand_descr, list_brand(nb)) & ind_yr;

    prod.p0 = mean(price(ind_nb));
    prod.x_rand = unique(x_rand(ind_nb, :), 'rows');
    prod.x_rand2 = unique(x_rand2(ind_nb, :), 'rows');
    prod.owner_ind = logical(unique(owner_ind(ind_nb, :), 'rows')); 
    prod.id = unique(brand_id2(ind_nb));
    prod.owner = unique(owner_name2(ind_nb));
    prod.brand = list_brand(nb);
    prod.ind_endogenous = true; % this product's price changes in the equilibrium
    prod.craft = unique(craftdummy(ind_nb));

    % placeholder options; not in use when we assume beer companies to set
    % retail prices
    prod.wp = 0; % wholesale prices
    prod.retail_ind = true; % retail index 
    

    % coordinate
    prod.brewer_id = unique(brewer_id(ind_nb));

    ind = brewer_loc_tbl.brewer_id == prod.brewer_id;
    prod.brewer_crdnt = [brewer_loc_tbl.latitude_brewer(ind), brewer_loc_tbl.longtitude_brewer(ind)];
    prod.brewer = unique(brewer_loc_tbl.brewer(ind));
    
    dist = nan(n_county_cty, 1);
    for nc = 1 : n_county_cty
        dist(nc, :) = min(lldistkm_vec(prod.brewer_crdnt, list_county_crdnt(nc,:)));        
    end
    prod.dist = dist;

    try
        if min(dist)<50 && ~ismember(prod.brewer, {'ab inbev', 'miller brewing co', 'coors brewing co'})
            prod.in_state = true;
        else
            prod.in_state = false;
        end
    catch err
        prod.brewer
        prod.brewer_id
        prod.owner
        prod.brand
        
        rethrow(err)
    end

    if ismember(prod.brewer, 'sierra nevada brewing co')
        prod.sierra = true;
    else
        prod.sierra = false;
    end

    if size(prod.x_rand, 1)~=1 || size(prod.owner_ind, 1)~=1 ...
            || size(prod.id, 1)~=1 || size(prod.owner, 1)~=1 || size(prod.brand, 1)~=1
        error('brand indexing error');
    end

    prod_attributes{nb} = prod;
end


%% compute distance*coeff, brand FE, mkt FE in both demand and mc and then draw them from their asymptotic distributions
list_geomkt = unique(geomkt);
ind_geomkt = dummyvar(geomkt); % n_obs x n_geo
n_geo = max(geomkt);

%(1) craft products
ind_s = craftdummy & ind_yr;
N_s = nnz(ind_s);

demand_resid_s = demand_resid(ind_s);
[list_brand_s, ~, id_brand_s] = unique(brand_descr(ind_s));

X_s = [dummyvar(id_brand_s), ind_geomkt(ind_s, :), dist_x(ind_s, :)]; 

%(A) regress (distance bin FE+brand FE+mkt FE) in demand on distance bin dummies, brand dummies, mkt dummies
b_demand_s = pinv(X_s'*X_s)*X_s'*demand_resid_s;
r_demand_s = demand_resid_s - X_s*b_demand_s;

qrs = quantile(r_demand_s, [0.01 0.99]);
r_demand_s(r_demand_s>qrs(2)) = qrs(2);
r_demand_s(r_demand_s<qrs(1)) =  qrs(1);

ind_s_sim = randsample(N_s, N_sim, true);
ind_s_mo_sim = randsample(12, N_sim, true);
r_demand_sim_s = r_demand_s(ind_s_sim) + para_mo(ind_s_mo_sim)';

%(B) draws from the asymptotic distributions
b_demand_s_boot = nan(N_boot, length(b_demand_s));
r_demand_sim_s_boot = nan(N_boot, N_sim);
for nb = 1 : N_boot
    % first stage statistical errors
    b_demand_s_nb = pinv(X_s'*X_s)*X_s'* demand_resid_boot(ind_s, nb);
    r_demand_s_nb = demand_resid_boot(ind_s, nb) - X_s*b_demand_s_nb;
    Vb_demand_nb = 1/N_s*X_s'*diag(r_demand_s_nb.^2)*X_s;

    % second stage statistical errors
    % sample from the asymptotic distribution
    b_demand_s_boot(nb, :) = mvnrnd(b_demand_s_nb, Vb_demand_nb);

    qrs_nb = quantile(r_demand_s_nb, [0.01 0.99]);
    r_demand_s_nb(r_demand_s_nb>qrs_nb(2)) = qrs_nb(2);
    r_demand_s_nb(r_demand_s_nb<qrs_nb(1)) =  qrs_nb(1);

    r_demand_sim_s_boot(nb, :) = r_demand_s_nb(ind_s_sim) ...
        + para_mo_boot(nb, ind_s_mo_sim)';
end

%(C) regress mc on  (distance bin FE), (brand FE), (mkt FE), (month FE)
mc_s = mc(ind_s);

ind_mc_mo_s = size(X_s, 2) + (1 : 11);
X_s = [X_s, mo_dummy(ind_s, 1:11)];

b_mc_s_all = pinv(X_s'*X_s)*X_s'*mc_s; % all FE including months
r_mc_s = mc_s - X_s*b_mc_s_all;

qm2s = quantile(r_mc_s, [0.01 0.99]);
r_mc_s(r_demand_s>qm2s(2)) = qm2s(2);
r_mc_s(r_demand_s<qm2s(1)) =  qm2s(1);

b_mc_s_mo = [b_mc_s_all(ind_mc_mo_s); 0]; % mont FE 
b_mc_s = b_mc_s_all(1:ind_mc_mo_s(1)-1); % all FE excluding month FE

r_mc_sim_s = r_mc_s(ind_s_sim) + b_mc_s_mo(ind_s_mo_sim);

%(D) draws from the asymptotic distributions
b_mc_s_boot = nan(N_boot, length(b_mc_s));
r_mc_sim_s_boot = nan(N_boot, N_sim);
for nb = 1 : N_boot
    % first stage statistical errors
    b_mc_s_all_nb = pinv(X_s'*X_s)*X_s'* mc_boot(ind_s, nb);
    r_mc_s_nb = mc_boot(ind_s, nb) - X_s*b_mc_s_all_nb;
    Vb_mc_nb = 1/N_s*X_s'*diag(r_mc_s_nb.^2)*X_s;

    % second stage statistical errors
    % sample from the asymptotic distribution
    b_mc_s_boot_draw = mvnrnd(b_mc_s_all_nb, Vb_mc_nb);
    b_mc_s_boot(nb, :) = b_mc_s_boot_draw(1:ind_mc_mo_s(1)-1);
    b_mc_s_mo_boot = [b_mc_s_boot_draw(ind_mc_mo_s)'; 0];

    qm2s_nb = quantile(r_mc_s_nb, [0.01 0.99]);
    r_mc_s_nb(r_mc_s_nb>qm2s_nb(2)) = qm2s_nb(2);
    r_mc_s_nb(r_mc_s_nb<qm2s_nb(1)) = qm2s_nb(1);

    r_mc_sim_s_boot(nb, :) = r_mc_s_nb(ind_s_sim) + b_mc_s_mo_boot(ind_s_mo_sim);
end

%(E) use the estimates and draws to construct draws of market-specific demand and mc shocks
for j = 1 : length(list_brand_s)
    j_list_id = find(ismember(list_brand, list_brand_s(j)));


    demand_resid_j = nan(n_geo, N_sim, N_boot+1);
    mc_j = nan(n_geo, N_sim, N_boot+1);

    brand_fe_j = zeros(1, length(list_brand_s));
    brand_fe_j(j) = 1;
    for m = 1 : n_geo
        mkt_fe_j = zeros(1, length(list_geomkt));
        mkt_fe_j(list_geomkt==m) = 1;

        ind = ismember(list_county, store_county_id(geomkt==m));
        dist_dummy_m = dist_dummy(prod_attributes{j_list_id}.dist(ind));

        % use nearest neighbor extrapolation by distance
        X_js = [brand_fe_j, mkt_fe_j, dist_dummy_m];

        demand_resid_j(m, :, 1) = X_js*b_demand_s + r_demand_sim_s;
        mc_j(m, :, 1) = X_js*b_mc_s + r_mc_sim_s;

        for nb = 1 : N_boot
            demand_resid_j(m, :, nb+1) = X_js*b_demand_s_boot(nb, :)'...
                + r_demand_sim_s_boot(nb, :);
            mc_j(m, :, nb+1) = X_js*b_mc_s_boot(nb, :)'...
                + r_mc_sim_s_boot(nb, :);
        end

    end
    prod_attributes{j_list_id}.demand_resid = demand_resid_j(:,:,1);
    prod_attributes{j_list_id}.demand_resid_boot = demand_resid_j(:,:,2:end);
    prod_attributes{j_list_id}.mc = mc_j(:,:,1); 
    prod_attributes{j_list_id}.mc_boot = mc_j(:,:,2:end); 
end

clear demand_resid_s list_brand_s id_brand_s X_s mc_s...
b_demand_s r_demand_s r_demand_sim_s_boot ...
b_demand_s_nb r_demand_s_nb Vb_demand_nb r_demand_s_nb  ...
b_mc_s_all r_mc_s b_mc_s_mo b_mc_s r_mc_sim_s r_mc_sim_s_boot ...
b_mc_s_all_nb r_mc_s_nb Vb_mc_nb ...
b_mc_s_boot_draw b_mc_s_mo_boot r_mc_s_nb ind_mc_mo_s ind_s_mo_sim ind_s_sim

%(2) other products
ind_o = ~craftdummy & ind_yr;
N_o = sum(ind_o);

demand_resid_o = demand_resid(ind_o);
[list_brand_o, ~, id_brand_o] = unique(brand_descr(ind_o));

X_o = [dummyvar(id_brand_o), ind_geomkt(ind_o, :), dist_x(ind_o,:)];

%(A) regress (distance bin FE+brand FE+mkt FE) in demand on distance bin dummies, brand dummies, mkt dummies
b_demand_o = pinv(X_o'*X_o)*X_o'*demand_resid_o;
r_demand_o = demand_resid_o - X_o*b_demand_o;

qro = quantile(r_demand_o, [0.01 0.99]);
r_demand_o(r_demand_o>qro(2)) = qro(2);
r_demand_o(r_demand_o<qro(1)) =  qro(1);

ind_o_sim = randsample(N_o, N_sim, true);
ind_o_mo_sim = randsample(12, N_sim, true);
r_demand_sim_o = r_demand_o(ind_o_sim) + para_mo(ind_o_mo_sim)';

%(B) draws from the asymptotic distributions
b_demand_o_boot = nan(N_boot, length(b_demand_o));
r_demand_sim_o_boot = nan(N_boot, N_sim);
for nb = 1 : N_boot
    % first stage statistical errors
    b_demand_o_nb = pinv(X_o'*X_o)*X_o'* demand_resid_boot(ind_o, nb);
    r_demand_o_nb = demand_resid_boot(ind_o, nb) - X_o*b_demand_o_nb;
    Vb_demand_nb = 1/N_o*X_o'*diag(sparse(r_demand_o_nb.^2))*X_o;

    % second stage statistical errors
    % sample from the asymptotic distribution
    b_demand_o_boot(nb, :) = mvnrnd(b_demand_o_nb, Vb_demand_nb);

    qro_nb = quantile(r_demand_o_nb, [0.01 0.99]);
    r_demand_o_nb(r_demand_o_nb>qro_nb(2)) = qro_nb(2);
    r_demand_o_nb(r_demand_o_nb<qro_nb(1)) = qro_nb(1);

    r_demand_sim_o_boot(nb, :) = r_demand_o_nb(ind_o_sim) ...
        + para_mo_boot(nb, ind_o_mo_sim)';
end

%(C) regress mc on  (distance bin FE), (brand FE), (mkt FE), (month FE)
mc_o = mc(ind_o);

ind_mc_mo_o = size(X_o, 2) + (1 : 11);
X_o = [X_o, mo_dummy(ind_o, 1:11)];

b_mc_o_all = pinv(X_o'*X_o)*X_o'*mc_o; % all FE including months
r_mc_o = mc_o - X_o*b_mc_o_all;

qm2o = quantile(r_mc_o, [0.01 0.99]);
r_mc_o(r_demand_o>qm2o(2)) = qm2o(2);
r_mc_o(r_demand_o<qm2o(1)) =  qm2o(1);

b_mc_o_mo = [b_mc_o_all(ind_mc_mo_o); 0];% month FE 
b_mc_o = b_mc_o_all(1:ind_mc_mo_o(1)-1);% all FE excluding month FE

r_mc_sim_o = r_mc_o(ind_o_sim) + b_mc_o_mo(ind_o_mo_sim);

%(D) draws from the asymptotic distributions
b_mc_o_boot = nan(N_boot, length(b_mc_o));
r_mc_sim_o_boot = nan(N_boot, N_sim);
for nb = 1 : N_boot
    % first stage statistical errors
    b_mc_o_mo_nb = pinv(X_o'*X_o)*X_o'* mc_boot(ind_o, nb);
    r_mc_o_nb = mc_boot(ind_o, nb) - X_o*b_mc_o_mo_nb;
    Vb_mc_nb = 1/N_s*X_o'*diag(sparse(r_mc_o_nb.^2))*X_o;

    % second stage statistical errors
    % sample from the asymptotic distribution
    b_mc_o_boot_draw = mvnrnd(b_mc_o_mo_nb, Vb_mc_nb);
    b_mc_o_boot(nb, :) = b_mc_o_boot_draw(1:ind_mc_mo_o(1)-1);
    b_mc_o_mo_boot = [b_mc_o_boot_draw(ind_mc_mo_o)'; 0];

    qm2o_nb = quantile(r_mc_o_nb, [0.01 0.99]);
    r_mc_o_nb(r_mc_o_nb>qm2o_nb(2)) = qm2o_nb(2);
    r_mc_o_nb(r_mc_o_nb<qm2o_nb(1)) = qm2o_nb(1);    
    r_mc_sim_o_boot(nb, :) = r_mc_o_nb(ind_o_sim) + b_mc_o_mo_boot(ind_o_mo_sim);
end

%(E) use the estimates and draws to construct draws of market-specific demand and mc shocks
for j = 1 : length(list_brand_o)
    j_list_id = find(ismember(list_brand, list_brand_o(j)));

    demand_resid_j = nan(n_geo, N_sim, N_boot+1);
    mc_j = nan(n_geo, N_sim, N_boot+1);

    brand_fe_j = zeros(1, length(list_brand_o));
    brand_fe_j(j) = 1;

    for m = 1 : n_geo
        mkt_fe_j = zeros(1, length(list_geomkt));
        mkt_fe_j(list_geomkt==m) = 1;
        ind = ismember(list_county, store_county_id(geomkt==m));
        dist_dummy_m = dist_dummy(prod_attributes{j_list_id}.dist(ind));

        % construct the demand and mc shocks
        X_jo = [brand_fe_j, mkt_fe_j, dist_dummy_m];
        demand_resid_j(m, :, 1) = X_jo*b_demand_o + r_demand_sim_o;
        mc_j(m, :, 1) = X_jo*b_mc_o + r_mc_sim_o;

        for nb = 1 : N_boot
            demand_resid_j(m, :, nb+1) = X_jo*b_demand_o_boot(nb, :)'...
                + r_demand_sim_o_boot(nb, :);
            mc_j(m, :, nb+1) = X_jo*b_mc_o_boot(nb, :)'...
                + r_mc_sim_o_boot(nb, :);
        end

    end
    prod_attributes{j_list_id}.demand_resid = demand_resid_j(:,:,1);
    prod_attributes{j_list_id}.demand_resid_boot = demand_resid_j(:,:,2:end);
    prod_attributes{j_list_id}.mc = mc_j(:,:,1);
    prod_attributes{j_list_id}.mc_boot = mc_j(:,:,2:end);
end
clear ind_o demand_resid_o list_brand_o id_brand_o X_o mc_o...
b_demand_o r_demand_o r_demand_sim_o_boot ...
b_demand_o_nb r_demand_o_nb Vb_demand_nb r_demand_o_nb  ...
b_mc_o_all r_mc_o b_mc_o_mo b_mc_o r_mc_sim_o r_mc_sim_o_boot ...
b_mc_o_all_nb r_mc_o_nb Vb_mc_nb ...
b_mc_o_boot_draw b_mc_o_mo_boot r_mc_o_nb ind_mc_mo_o ind_o_mo_sim ind_o_sim

%% simulate min and max change in variable profit profits
no_entry_mkt = find(sum(ind_geomkt(ind_s, :))<25); clear ind_s; %markets with fewer than 25 craft/month combinations are excluded
craft_mkt = setdiff(geomkt(ind_yr), no_entry_mkt);

N_mkt = length(craft_mkt);
N_dev = nnz(~is_fixed_prod);

endo_dev_ind = find(~is_fixed_prod);

prod_entry_ind = false(N_mkt, length(list_brand)); % the entry conditions of all products
entry_sim = nan(N_mkt, N_dev);
brand_id_sim = nan(N_mkt, N_dev);
owner_id_sim = nan(N_mkt, N_dev);
dist_sim = nan(N_mkt, N_dev); 

% when all other products are in the market
dev_profit_all = nan(N_mkt, N_dev);

% when none of the other products is in the market
dev_profit_one = nan(N_mkt, N_dev);

% simulated profits taking into account statistical errors
dev_profit_all_boot = nan(N_mkt, N_dev, N_boot);

dev_profit_one_boot = nan(N_mkt, N_dev, N_boot);

inc_profit = nan(N_mkt, 1);
pop_profit = nan(N_mkt, 1);


parfor m = 1 : N_mkt
    disp(['m=' num2str(m)]);

    % create a product attribute matrix specific to the market
    prod_attributes_m = prod_attributes;

    % choose the demand residuals and mc shocks specific to the marke
    for j = 1 : length(prod_attributes)
        prod_attributes_m{j}.demand_resid = prod_attributes{j}.demand_resid(craft_mkt(m), :);
        prod_attributes_m{j}.mc = prod_attributes{j}.mc(craft_mkt(m), :);
    end

    % identify the market-year
    ind_m = (geomkt==craft_mkt(m) & ind_yr);

    % find the brands in the market that are present for more than cutoff_month_count of months
    brand_tab_m = tabulate(brand_descr(ind_m));    
    count_brands_m = cell2mat(brand_tab_m(:, 2)); % count the number of months of a product in the year in ind_yr    
    prod_entry_ind_m = ismember(list_brand, brand_tab_m(count_brands_m>cutoff_month_count, 1));% consider these brands as being in the market
    prod_entry_ind(m, :) = prod_entry_ind_m; %entry decisions of all products
    
    entry_sim(m, :) = prod_entry_ind_m(endo_dev_ind); %entry decisions of non-fixed products

    % market level characteristics
    inc_profit(m) = unique(inc_mkt(ind_m));
    pop_profit(m) = mean(unique(mktsize(ind_m)));

    % market specific price and income distribution
    sim_pcoeff_t = unique(price_coeff(ind_m,:), 'rows');
    loginc_t = unique(loginc(ind_m, :), 'rows');

    for nd = 1 : N_dev
        % brand id, owner id, and distance
        brand_id_sim(m, nd) = prod_attributes_m{endo_dev_ind(nd)}.id;
        owner_id_sim(m, nd) = find(prod_attributes_m{endo_dev_ind(nd)}.owner_ind);
        ind = ismember(list_county, store_county_id(ind_m));
        dist_sim(m, nd) = prod_attributes_m{endo_dev_ind(nd)}.dist(ind);

        % product configuration where all other endogenous products are in
        prod_entry_ind_m_all = prod_entry_ind_m;
        prod_entry_ind_m_all(~is_fixed_prod) = true; % this product is in (i.e., j=1, -j=1)

        prod_entry_ind_m_all_base = prod_entry_ind_m_all;
        prod_entry_ind_m_all_base(endo_dev_ind(nd)) = false; %this product is not in (i.e., j=0, -j=1)

        % product configuration where none of the other endogenous products are in
        prod_entry_ind_m_one = prod_entry_ind_m;
        prod_entry_ind_m_one(~is_fixed_prod) = false;

        prod_entry_ind_m_one_base = prod_entry_ind_m_one; % this product is not in  (i.e., j=1, -j=0)
        prod_entry_ind_m_one(endo_dev_ind(nd)) = true; % this product is in (j=1, -j=0)

        % each firm's profit at (j=0, -j = 1)
        oem_profit_all_base = equi_ind(prod_entry_ind_m_all_base, prod_attributes_m, [], sim_pcoeff_t, ...
            setup, para_nonlinear, randomdraw, ...
            loginc_t, logincnorm);

        % each firm's profit at (j=1,-j=1)
        oem_profit_all_add = equi_ind(prod_entry_ind_m_all, prod_attributes_m, [], sim_pcoeff_t, ...
            setup, para_nonlinear, randomdraw, ...
            loginc_t, logincnorm);

        % each firm's profit at (j=0, -j=0)
        if any(prod_entry_ind_m_one_base)
            oem_profit_one_base = equi_ind(prod_entry_ind_m_one_base, prod_attributes_m, [], sim_pcoeff_t, ...
                setup, para_nonlinear, randomdraw, ...
                loginc_t, logincnorm);
        else
            oem_profit_one_base = zeros(size(oem_profit_all_add));
        end

        % each firm's profit at (j=1, -j=0)
        oem_profit_one_add = equi_ind(prod_entry_ind_m_one, prod_attributes_m, [], sim_pcoeff_t, ...
            setup, para_nonlinear, randomdraw, ...
            loginc_t, logincnorm);

        % change in profit of j's firm when -j=1
        dev_profit_all(m, nd) = (oem_profit_all_add(prod_attributes_m{endo_dev_ind(nd)}.owner_ind) - ...
            oem_profit_all_base(prod_attributes_m{endo_dev_ind(nd)}.owner_ind))*pop_profit(m);

        % change in profit of j's firm when -j=0
        dev_profit_one(m, nd) = (oem_profit_one_add(prod_attributes_m{endo_dev_ind(nd)}.owner_ind) - ...
            oem_profit_one_base(prod_attributes_m{endo_dev_ind(nd)}.owner_ind))*pop_profit(m);



             
        % repeat the above computation for each bootstrap draws
        for nb = 1 : N_boot
            prod_attributes_m_nb = prod_attributes;

            for j = 1 : length(prod_attributes)
                prod_attributes_m_nb{j}.demand_resid = prod_attributes{j}.demand_resid_boot(craft_mkt(m), :, nb);
                prod_attributes_m_nb{j}.mc = prod_attributes{j}.mc_boot(craft_mkt(m), :, nb);
            end

            oem_profit_all_base = equi_ind(prod_entry_ind_m_all_base, prod_attributes_m_nb, [], sim_pcoeff_t, ...
                setup, para_nonlinear, randomdraw, ...
                loginc_t, logincnorm);

            oem_profit_all_add = equi_ind(prod_entry_ind_m_all, prod_attributes_m_nb, [], sim_pcoeff_t, ...
                setup, para_nonlinear, randomdraw, ...
                loginc_t, logincnorm);

            if any(prod_entry_ind_m_one_base)
                oem_profit_one_base = equi_ind(prod_entry_ind_m_one_base, prod_attributes_m_nb, [], sim_pcoeff_t, ...
                    setup, para_nonlinear, randomdraw, ...
                    loginc_t, logincnorm);
            else
                oem_profit_one_base = zeros(size(oem_profit_all_add));
            end

            oem_profit_one_add = equi_ind(prod_entry_ind_m_one, prod_attributes_m_nb, [], sim_pcoeff_t, ...
                setup, para_nonlinear, randomdraw, ...
                loginc_t, logincnorm);

            dev_profit_all_boot(m, nd, nb) = (oem_profit_all_add(prod_attributes_m_nb{endo_dev_ind(nd)}.owner_ind) - ...
                oem_profit_all_base(prod_attributes_m_nb{endo_dev_ind(nd)}.owner_ind))*pop_profit(m);

            dev_profit_one_boot(m, nd, nb) = (oem_profit_one_add(prod_attributes_m_nb{endo_dev_ind(nd)}.owner_ind) - ...
                oem_profit_one_base(prod_attributes_m_nb{endo_dev_ind(nd)}.owner_ind))*pop_profit(m);

            all_prod_profit_all_boot(m, nd, nb) = ...
                oem_profit_all_add(prod_attributes_m_nb{endo_dev_ind(nd)}.owner_ind)*pop_profit(m);

        end
    end
end


%% covariates in FC functions
% in state and non-ABI/MillerCoors breweries
in_state_sim = zeros(N_mkt, N_dev);
craft_sim = zeros(N_mkt, N_dev);


list_brand_id = nan(length(list_brand), 1);

for j = 1 : length(list_brand)
    list_brand_id(j) = unique(brand_id2(ismember(brand_descr, list_brand(j))));
end

for m = 1 : N_mkt
    for nd = 1 : N_dev
        prod_m_nd = prod_attributes{list_brand_id == brand_id_sim(m, nd)};

        % goose island has been bought
        if isequal(prod_m_nd.brewer, 'goose island beer co')
            craft_sim(m, nd) = 0;
        else
            craft_sim(m, nd) = prod_m_nd.craft;
        end
        in_state_sim(m, nd) = prod_m_nd.in_state;

    end
end

%% save
% will be used later
varlist_for_next_step ={'prod_attributes','dev_profit_one','dev_profit_all','dev_profit_one_boot','dev_profit_all_boot', ...
    'in_state_sim','craft_sim','dist_sim','entry_sim','pop_profit','endo_dev_ind', ...
    'para_nonlinear','randomdraw','logincnorm',...
    'yr','craftdummy','quantity','geomkt','store_county_id','loginc',...
    'owner_name2','fixed_brand','fixed_craft_brands','list_brand','brewer_id','brewer_name','brand_descr',...
    'craft_mkt','list_county','price_coeff','setup','cutoff_month_count', 'N_boot'};
save([output_folder, 'fc_covariates', num2str(supply_est_yr)], varlist_for_next_step{:}, '-v7.3');